%%
clear; clc;

omegaX = 2*pi* 0.0; omegaY = 2*pi* 10.0; omegaZ = 2*pi*34e3;
kc = 1.0; omegaC = kc *omegaZ; Rmax = 10e3; T1=inf; T2=10e-6;

gamma = Rmax + 1/T1 + 1/T2;

x = 3;

%%
bx=omegaX/gamma; by=omegaY/gamma; z0=omegaZ/gamma; zc=omegaC/gamma; 

t=linspace(0, 2*pi/zc, 100); 

% analytic solution
import Algorithm.BlochEquation.ParametricBlochEquation

m = ParametricBlochEquation(T1, T2, ...
                           'ratio_x', x, 'ratio_kc', kc);

m.set_Omega('omegaX', 2*pi* 0.0,  ... 
            'omegaY', 2*pi* 10.0, ...
            'omegaZ', 2*pi* 34e3) ...
 .set_pumping( Rmax )...
 .solveAnalytically();

sol = m.getSolution;


% numeric solution
[mx_num, my_num, model] = num_sol(T1, T2, omegaX, omegaY, omegaZ, Rmax, omegaC, x, t);

% plot
subplot(2, 1, 1)
plot(t, sol.mx, 'rd', t, sol.my, 'bo', t, mx_num, 'r-', t, my_num, 'b-');
legend({'mx', 'my'});
xlim(minmax(t))
subplot(2, 1, 2)
plot(t, (sol.mx-mx_num)./mx_num, 'rd-', t, (sol.my-my_num)./my_num, 'bo-')
xlim(minmax(t))

%%
function [mx, my, model] = num_sol(T1, T2, omegaX, omegaY, omegaZ, Rmax, omegaC, x, t)
    tmax = 100.0/(Rmax + 1/T1 + 1/T2);

    model = Algorithm.BlochEquation.BlochEquation(T1, T2, ...
                                    'steadyOpt', 'manual', ...
                                    'tRange', [0 tmax], ...
                                    'steady',    [0 0 -Rmax].');

    model.set_Omega('omegaZ', -model.cosine_pump(x * omegaC, omegaC/2/pi, omegaZ),...
                'omegaX', omegaX, ...
                'omegaY', omegaY ).set_pumping( Rmax ).solve();


%     model.plot_solution;
    sol=model.getSolution;

    period = 2*pi/omegaC;
    nP = floor(tmax/period);
    tlist = linspace( (nP-1)*period, nP*period, length(t) );

    mx = sol.mx(tlist);
    my = sol.my(tlist);

end